Constructing the generalized Gibbs ensemble after a quantum quench 
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Using a numerical renormalization group based on exploiting an underlying exactly solvable non- 
relativistic theory, we study the out-of-equilibrium dynamics of a ID Bose gas (as described by the 
Lieb-Liniger model) released from a parabolic trap. Our method allows us to track the post-quench 
dynamics of the gas all the way to infinite time. We also exhibit a general construction, applicable 
to all integrable models, of the thermodynamic ensemble that has been suggested to govern this 
dynamics, the generalized Gibbs ensemble. We compare the predictions of equilibration from this 
ensemble against the long time dynamics observed using our method. 
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Understanding non-equilibrium quantum quench be- 
havior in low-dimensional systems is a difficult theoret- 
ical challenge. Because one is initializing the system in 
a state that is not an eigenstate, this behavior is deter- 
mined not merely by the system's ground state (or a small 
number of excited states), but rather by some coherent 
sum of a large number of eigenstates. If one wants to ex- 
plore the emergence of a resulting steady state, the time 
evolution of this coherent sum must then be tracked over 
long periods of time. This problem confronts theorists 
who wish to understand dynamics in perturbed quantum 
gases [li i2j , ultrafast phenomena in superconductors i3j , 
and questions of thermalization in integrable systems [4] . 

This last set of questions arise because of the surprising 
experimental finding that a perturbed one-dimensional 
Bose gas retains memory of its initial non-equilibrium 
state over long periods of time [1 and does not appear 
to relax to a state of thermodynamic equilibrium. To un- 
derstand this, it was proposed [1] that equilibriation does 
occur but not as described by a grand canonical ensem- 
ble (GCE). Instead the ensemble describing equilibria- 
tion needs to take into account the additional, non-trivial 
conserved quantities that, at least according to the the- 
oretical minimal model of the gas (the Lieb-Liniger (LL) 
model [5]), are present in the system. This new ensemble 
has been dubbed the generalized Gibbs ensemble (GGE). 
The GGE takes as the density matrix 



PGGE 



= Z ^ exp(- 



/^ 



P^Q^) 



(1) 



where the Qi form an independent, complete sequence of 
conserved quantities in the system and Pi correspond to a 
set of generalized (inverse) temperatures. Computation 
of this density matrix is non-trivial and has only been 
successfully accomplished in certain special limits. Most 
of these limits are in models where interactions (though 
not necessarily correlation functions) correspond to a free 
model (the hard core limit of the interacting Bose gas 
[3], quadratic Hamiltonians f^, Luttinger liquids fTj, the 



sine-Gordon model at the free-fermion point and in the 
semi-classical limit [8 , and the quantum Ising model in 
the absence of a longitudinal field [21 [TU]). A notable 
exception was the study of Fioretto and Mussardo [TT] 
where it was possible to study quenches in general in- 
teracting integrable models but with the restriction to a 
very special set of quench protocols. 

It is against this backdrop that we present a general 
methodology able to study non-equilibrium behavior and 
quench dynamics of low-dimensional interacting models, 
both integrable and non-integrable. This method is pred- 
icated on a numerical renormalization group (NRG) able 
to study models which can be represented as perturbed 
integrable and conformal field theories (GET) [I^: 



H = H^ 



Integrable /C FT 



V, 



perturbation ■ 



(2) 



The LL model in a trapping potential takes this form. 
We believe that this methodology is a valuable addition 
to other general methodologies used to study dynamics 
in low-dimensional systems such as the time-dependent 
density matrix renormalization group |13H17j . At least 
for a subset of quenches, where we quench into an inte- 
grable system (say by turning off the trapping potential 
in a LL system), we can track the dynamics for all times. 

Goncomitant with the introduction of this tool to study 
quench dynamics, we present a general methodology to 
compute the density matrix of the GGE using informa- 
tion arising from the application of the NRG. We show 
how one can write down a simple set of equations govern- 
ing the GGE and how the entire infinite set of generalized 
temperatures, {/3i}^i can be readily determined. 

The specific example we consider is the LL model per- 
turbed by a one-body parabolic trap V(x) = muj'^x^ 12, 



N 



H = - 



2m ^-^ dx^ 



2c^<5(a 



x,) + Y.Vix,), (3) 



(we will work in units where 2m = h = 1). In running 
the NRG, we use the basis of eigenstates of the LL model 



and their matrix elements with respect to the trapping 
potential. Both the description of the states and the com- 
putation of matrix elements in the LL model are much 
more complicated than the examples of relativistic field 
theories where the NRG has been applied previously. The 
states in the LL model consist of N strongly interacting 
particles and not few-particle excitations above the true 
vacuum state, while the matrix elements do not see a 
chiral factorization as in a relativistic gapless theory but 
are N-dimensional determinants [TH]. To tackle this, we 
took recourse to a highly optimized set of routines known 
as ABACUS ^^ which solves and evaluates all equations 
needed to characterize both the necessary eigenstates and 
their matrix elements. This package has been shown to 
be able to successfully compute dynamical response func- 
tions for the LL model 19 . 

We first use the NRG to extract the ground state of the 
LL model in a trap ^0\ . The NRG produces the ground 
state of the gas, \iP)gSj ^^ a linear combination of exact 
eigenstates, \s), of the LL model: |V')gs = Ss^sls). In 
order to accurately describe the ground state in the NRG 
procedure we typically consider on the order of 10^ — 
10^ states. We then consider a sudden release of the 
trap, that is we will study the gas where we quench into 
an integrable model. For these types of quenches our 
methodology gives us the ability to study the evolution 
of the gas for arbitrary times. Each state, \s), appearing 
in the ground state is characterized by a set of N (one 
for each particle) rapidities (quasi-momenta) {Xn}n=i- 
These rapidities are solutions to the Bethe equations, 
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(4) 



and can be readily obtained to arbitrary accuracy. With 
the NRG we can compute the coefficients c^ with rea- 
sonably high accuracy [20'. Time evolution under the 
post-quench Hamiltonian (the unperturbed LL model) is 
extremely simple, li Eg is the energy of state |s), the time 
evolution is described by \tp{t))GS = Ss Cse~*'^=*|s). Be- 
cause each state's energy, E^, is given in terms of the A„'s 
as ^^ A^, we can compute the phases appearing in the 
above sum to arbitrary accuracy for arbitrary time. 

To characterize the evolution of the gas in the long time 
limit we compute the momentum distribution function 
(MDF) rife = (V'feV'fc) in the diagonal ensemble (DE). An 
observable O^O in this ensemble is simply given by 



(O^O) 



DE 



El 



(s\0^0\s) 



(5) 



To compute this correlation function we insert a resolu- 
tion of the identity between O^ and O and use a specially- 
designed version of ABACUS for excited states to com- 
pute all of the necessary matrix elements [5^ . 

In Fig. [T] we plot the MDF in the DE of the gas post- 
release for two values of c (c = 10 and c = 7200) and 
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FIG. 1: The MDF in the DE of the gas after release from a 
trap for c = 10 (top) and c = 7200 (bottom). Shown are the 
gases at (iV = L = 14,0; = 0.64), (iV = L = 28,aj = 0.32), 
and (A'' — L — 56, to — 0.16). Error bars are given for the 
A'' = L = 56 data alone and are estimated from the speed 
of convergence of the NRG (see [201) ~ we believe the A'^ = 
L = 14, 28 data is completely converged. The MDF of the 
untrapped gas (A'^ = L = 56) is shown for comparison as is 
the analytic expression available for the Tonks-Girardeau gas 
c = oo from Ref. '21 . 



for a variety of system sizes, with uL fixed and keeping 
N — L. For comparison we also plot the MDF of the gas 
in its ground state. 

We see, as expected, that the MDF of the gas is per- 
turbed from that of the ground state at low momenta but 
remains unchanged from the ground state MDF at higher 
momenta. The relative insensitivity to different values of 
iV = L, w is consistent with a perturbative (in uj) compu- 
tation of the MDF in the DE at c = oo which shows 
n{k)DE = nik)Gs + i&T'-tffl^ + ^(^')- Here 
n{k)Gs is the MDF of the ground state, the constant 
Bq « 0.5124 [22j, and vp is the velocity of the gas. The 
scaling with A^, L, and oj indicated by this expression im- 
plies that variations in n{k)£)E between different system 
sizes in Fig. 1 are due to finite size corrections which are 
small (on the order of the symbol size) . As an important 
check of our results, the high momenta tails of the MDF's 
at c = 7200 behave as the predicted fc""* pTlEHl l^. 



While the diagonal ensemble tells us what the final 
steady state of the gas is after its release, a question 
of primary interest is whether the steady state can be 
associated with some ensemble. It has been postulated 
[3] that for a quench into an integrable system the correct 
ensemble to use is the GGE ensemble in Eqn. [l] The Qi's 
are here non-trivial polynomials in the field operators 
(and their derivatives) [25]. The action of the QiS on 
the states, \s) is straightforward. With each state, |s), 
characterized by a set of A^ rapidities, A^, the action of the 
Qi upon \s) is Qi\s;\i,-- ■ ,Xn) = Sj A* |s; Ai, • • • ,Xn), 
that is to say, Qi acts on the state like an i-th power 
sum. This shows that the Qi's are both a complete and 
independent set of charges inasmuch as the polynomials 
form a complete and independent basis in the space of 
single variable functions. 

To compute pgge the most straightforward path is 
to compute (Qi) at t — and insist that the set 
of /3i's is such that TripcGEQi) gives the same an- 
swer. In the case of the hard core limit this is read- 
ily doable as the Qi's can be written in terms of a 
more amenable basis, the momentum occupation num- 
bers: Qi — X^A'^*'^^' where n\ tells you whether there 
is a particle with rapidity of the form A — 2T:m/L for 
m e Z. In this basis of charges, {uxjcge simplifies to 
Tr(exp(— /3AnA)'T-A)/Trexp(— /3AnA), i-O- for such expec- 
tation values the ensemble factorizes, and Px is readily 
computed. This simplification, however, does not ex- 
ist away from the hard core limit and we are instead 
left with a complicated non-linear minimization problem 
which on the face of it does not obviously have a solu- 
tion. We now show that it does and that the PiS can 
be computed readily. We do so through a (generalized) 
thermodynamic Bethe ansatz |26| . 

Because the action of the charges Qi on the states, |s), 
are given simply in terms of the rapidities, A^ , identifying 
the state, to ask that {Qi)t=o — {Qi)GGE amounts to 
asking whether there is a set of A's, {Aj}^ 



■^]\jLi, such that 



{Q^ 



Ea^ 



1,2,' 



There is in fact such a set. We can moreover determine its 
rapidity distribution, which we will call PggeW, directly 
from \iP)gs- To each state, \s; A^i, • • • , Xsn), we associate 
a distribution, Ps(A), governing the A's of that particular 
state: ps{X) = xJ2iH^- ^sz)- Then pgge{X) is the 
weighted sum of the ps(A)'s: 



Pgge{X) = ^1 



%(A). 



In particular ^ dXpGGE{X)X^ = L{Q^)t=o. 

Pgge contains, implicitly, all the information to char- 
acterize the action of pgge on a eigenstate of the LL 
model [2^. A distribution of A's must be consistent with 
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FIG. 2: eo(A) and p{\) for both the GGE and GGE ensembles 
for a gas with N = L — 56, c = 7200, and a prequench trap 
strength, oj = 0.256. For the GGE ensemble the effective tem- 
perature is r = 1.54. The quantities plotted are symmetric 
about A = 0. 



the Bethe equations (Eqn. HI). In the continuum limit, 
these equations can be rewritten as O [57] 



Pgge{X) + p'^ceW ^ 7r + 



dX' 

2^ 



K{X-X')pgge{X), 

(6) 

where PggeW ^^ the density of holes in the A- 
distribution and K{X) = 2c/{c^ + A^). Now the GGE 
is derived by the same principles as the grand canoni- 
cal ensemble: namely entropy is maximized subject to 
the constraints of fixed conserved charges (energy for 
the grand canonical ensemble, all the charges, Qi, for 
the GGE). Thus associated with GGE is a generalized 
free energy Fgge = JdXpGGE{X)eo^GGE{X) - S, where 
£o-gge{X) = J2il^i^^ i^ ^ generalized energy. It corre- 
sponds to the action of pgge on a state |s; Ai, • • • , Xn): 



Pgge\s; Ai, 



,A 



= — J2i ^O-GGsi^i) 



N 



-\s; Ai, • • • , Xn)- 



(7) 
In particular knowing Eq-gge then allows us to compute 
general expectation values in the GGE. While Eq-gge 
differs from its form in the grand canonical ensemble, 
S is the standard entropy [27] of a system with a given 
distribution of particles, Pgge^ a-nd holes, Pqq^- 



S 



dX 



{PGGE + Pgge) log (pggb + Pgge) 



-PGGE log PGGE - Pgge log Pgge 



(8) 



We now show that we can express Eq-gge in terms of 
Pgge that we derived from \iP)gs- 

If we minimize the generalized free energy we arrive at 
a constraint between the particle and hole distributions 
and Eo-GGE- 



e{X)=eo-gge{X)- 



dX' 

2^ 



X(A-A')log(l + e-^(^)), (9) 
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FIG. 3: The MDF (for iV = L = 56) in the GCE and GGE 
for the gas after release from a trap of strength a; = 0.16 for 
c = 10 (top) and c = 7200 (bottom). We again show the 
MDF of the untrapped gas for comparison (blue stars). 



where e = \o%(^p%Q^I pgge)- Thus to (ietermine Eq-GGB 
we take our knowledge of pggeW obtained from \'iP)gSi 
use Eqn. (pi) to determine Pqq^ which then gives us £(A). 
From Eqn. (lo]), we then can fix Eq-gge- 

Following this procedure we plot in Fig. [2] pgge and 
Sq^gge for the gas in the hard core limit. For comparison 
we plot what these quantities would be if instead of a gen- 
eralized Gibbs ensemble, the thermodynamics was gov- 
erned by the grand canonical ensemble. (In this case we 
use the standard thermodynamic Bethe ansatz equations 
[27] to determine what pgce and Eq-gce = /3(A^ — /i) 
need to be, i.e. what the effective temperature needs to 
be, if they are to reproduce the correct density and av- 
erage energy of the system Gs(V'l-f^lV')GS-) We see that 
both Pgge and Eo-ggb have considerably more struc- 
ture than that of their grand canonical counterparts. 

We now use this ability to compute £o-ggb(A), to 
compute various expectation values of observables in the 
GGE. In Fig. [3]we plot the MDF as computed in the DE 
and in both the GGE and GCE. The error estimate is 
computed similarly as in Fig. fl] (see [20] for details) . For 
the data at hand, we see that for low momenta the two 
ensemble averages, GGE and GCE, disagree with the DE. 



However the GGE provides a considerably better match 
to the DE than does the ordinary thermal ensemble GCE. 
From the finite size comparison (see Fig. 3 of [20], it can 
be argued (although not conclusively) that at small but 
finite fc, this difference will vanish with increasing system 
size. 



The disagreement between ensembles in the data is 
not entirely surprising. The logic of the GGE is such 
that it is expected to describe correlations that are local 
in space (and that involve a distance scale significantly 
smaller than the system size) . We thus do not expect the 
correlations at /c ~ 1/L to be particularly well described 
by the GGE. However there is the possibility that the 
differences between ensembles will remain at finite k > 
1/L even in the infinite volume limit. In recent work 
[25] the entropy associated with the DE was shown to 
be considerably smaller than that of the GGE implying 
that the DE is more tightly constrained than the GGE, 
i.e. the GGE seems to be missing correlations. It would 
be interesting to understand if this missing entropy is 
solely associated with non-local correlations. 



In conclusion, we demonstrated how an NRG based 
on exploiting the integrability of the LL model can be 
used to study the time-dependent evolution after a quan- 
tum quench where a ID gas is released from a parabolic 
trap. We have also demonstrated how to use the informa- 
tion arising from the NRG to construct the corresponding 
GGE which has been suggested as a possibility for gov- 
erning the post-quench dynamics. While we have focused 
on the LL model, this methodology is applicable to any 
non-relativistic integrable theory of which the Heisenberg 
and XXZ spin chains are two prominent examples. 
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Supplementary Material 



Description of the numerical renormalization group 



The NRG we employ is one appropriate to the study 
of continuum field theories P^ [5^ which in turn is based 
upon Wilson's numerical renormalization group first used 
to study quantum impurity problems [30] • The basic idea 
behind Wilson's numerical approach is that one performs 
a series of numerical diagonalizations which are ordered 
such that "important" states in the considered Hilbert 
space of the problem are taken into account first, while 
the effect of less important states are included only in 
subsequent diagonalizations. The sequence of diagonal- 
izations, the sequence of renormalizations, are done such 
that the numerical burden is the same at each step in the 
sequence. The metric that determines the order in which 
states are taken into account, however, is arbitrary and 
can be chosen to be appropriate for the problem at hand. 
In Wilson's case, namely the Kondo problem where the 
impurity spin sits at the end of a half line lattice, states 
involving only the Kondo spin together with nearby lat- 
tice sites are taken into account first. That such states 
are most significant is guaranteed by lattice hopping pa- 
rameters that decrease in magnitude the further one gets 
away from the impurity. In the case of a quantum crit- 
ical Ising model perturbed by a magnetic field (consid- 
ered in [12]) the Hilbert space used to form the matri- 
ces in the sequence of numerical diagonalizations is that 
of the quantum critical Ising model. The states in this 
model are ordered in terms of energy (relative to the un- 
perturbed theory). Here low-energy states are the most 
important as the spin operator coupling to the magnetic 
field is highly relevant and so are taken into account first 
by the renormalization group. 

However for models like the LL model in a trapping po- 
tential the complexity of the eigenstates (as we have dis- 
cussed in the main body of the text) means that energy 
alone is not a sufficient metric for the NRG to distin- 
guish important from less important states; using such 
a limited metric would make the procedure drastically 
sub-optimal. To overcome this problem we introduce a 
variational metric in the space of states similar to that 
used to compute the single particle spectrum of semi- 
conducting carbon nanotubes [29]. Although this latter 
problem could be represented as a perturbation of a con- 
formal field theory (CFT), the CFT was complex enough 
(four bosons) that the same issues arose. This varia- 
tional metric uses an iterative process which amounts 
to performing successively higher order computations in 
perturbation theory to determine which states are likely 
to significantly contribute to the low-energy spectrum of 
the fully perturbed theory. 
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FIG. 4: Comparing analytics with NRG-fABACUS for 
computation of the ground state energy for different trap 
strengths and values of c at N = L = 56. Left: Ground 
state energy as a function of oj. Right: Ground state energy 
at different values of c - analytics include up to 1/c^ contri- 
butions. Numerical data for different u have been offset from 
each other. 



Description of the gas in the trap in equilibrium 

While the primary concern of the accompanying letter 
was the discussion of correlations, we feel it is important 
here to provide some details of the method so that readers 
can be reassured that we can describe the gas in the trap 
accurately, that is to say, that we can provide a reason- 
able description of the pre-quench state. We will provide 
a more detailed description of these computations in a 
future publication |31j . 

To demonstrate this we focus on the large c limit where 
in the extreme Tonks-Girardeau limit (c ~ oo) the system 
reduces to one of free fermions (at least for the compu- 
tation of the energy) in a trap and where 1/c corrections 
can be systematically computed by mapping the Bose 
system to one of fermions interacting with an ultra-short 
ranged potential with strength proportional to 1/c (33) . 

We present the computation of the ground state ener- 
gies in Fig. 4. In the left panel we present the ground 
state energies of the gas (with N = L = 56) in traps 
of different strengths, w. We get good agreement be- 
tween the NRG computation and the analytics (better 
than 0.02% for the first three trap values and about 1% 
for the largest of the trap values, a; = 0.256, studied). In 
the right panel we show the computation of the ground 
state energies as a function of c for values running from 
c = 7200 to c = 10 for the same four values of the trap 
strength. In order to match analytics with numerics we 
needed to include corrections up to 1/c'^. Computing the 
1/c correction is straightforward. The 1/c^ correction, 
when computed naively with second order perturbation 
theory, shows an ultraviolet divergence related to the 
ultra short ranged potential of the equivalent fermionic 
model. This divergence can however be regulated with a 
point splitting procedure adopted from Ref. |34j . How- 
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FIG. 5; Density profiles in the trapped gas at two different 
values of c comparing NRG+ ABACUS with analytics for N = 
L = 56: top) uj = 0.16; bottom) uj = 0.183. 



ever these first two corrections are not enough to get 
good agreement for the c = 10 values of the ground state 
energies. We thus include the 1/c^ correction coming 
from the untrapped gas and which is readily computed 
from its integrabihty [5] . With these three analytic con- 
tributions in place, we see we get excellent agreement for 
the first three trap values while the strongest trap value 
(w = 0.256) continues to see a deviation on the order of 
1% for all values of c. 

Finally we consider our ability to accurately compute 
the density profile of the gas in the trap. This is a much 
more complicated quantity to compute than the ground 
state energies, since the matrix elements of the density 
operator must be employed. In Fig. 5 we plot the density 
profile of the gas for two values of trapping strength and 
two values of c. Again we work at larger values of c where 
we can compare our numerics to an analytic computation 
(c = oo plus first order 1/c corrections). We see that we 
get good agreement in all cases. 



Finite Size Corrections to the MDF 

In Fig. 6 we present data for system sizes N = L = 28 
and iV = L = 56 of the differences of the MDF between 
the GGE and the DE as well as the GCE and the DE. 
We see that generically differences are present for small 
momenta. However these differences seem to behave dif- 
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FIG. 6: The difference between the MDF as computed in the 
DE and as computed in the GGE and GCE for A^ = L = 28 
(left) and A'^ = L = 56 (right). Top: The gas after release 
from a trap of strength uj — 0.32 (TV — L — 28) and lj = 
0.16 {N = L = 56) for c = 7200. Bottom: The gas after 
release from the same traps but for c = 10. The GGE estimate 
is seen to be more accurate than the GCE throughout the 
range of values considered. 



ferently for the different ensembles. For DE-GCE these 
differences remain approximately the same as one moves 
from N = L ^ 28 to N = L ^ 56. However for DE- 
GGE, these differences are notably reduced in increasing 
the system size, particularly for the c = 10 case. While 
we do not have sufficient data to perform a finite size 
scaling analysis, it appears that the difference DE-GGE 
is vanishing as system size grows, while the difference 
DE-GCE remains at some finite, though k-dependent, 
value. 



The SSF in the Various Ensembles 

While we focused on the MDF in the main body of the 
text, we also have computed the density-density corre- 
lation function (static structure factor (SSF)) SPP{k) = 



iip^ 



kP-k) (with pk = J2, 



^4+c 



V'9 



- fc+o -r H for definitions of the 
density operator see Ref. [32)- In Fig. [t] we plot the 
SSF in the DE for c = 10 and c = 7200. We see that the 
effects of the trap are restricted to small momenta and 
that they are enhanced as one reduces c, moving away 
from the hardcore limit. 

The role of scaling with system size is more compli- 
cated with the SSF than with the MDF. With the MDF 
we were able to argue (at least for weak values of the 
trap) that it remained invariant under the following scal- 
ing: N,L,cu ^ 2N, 2L, w/2. This is not true for the SSF. 
At weak values of to, we find that the SSF in the DE is 
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FIG. 7: The SSF in the DE of the gas after release from a 
trap of strength uj for c = 10 (top) and c = 7200 (bottom). 
Shown are the gases at (A'' = L = 14,tj = 0.64), {N ^ L = 
28, uj^ 0.32), and (N = L = 56,lo = 0.16). The N ^ L ^ 56 
data for higher momenta k ~ 2fcF is not fuUy saturated. The 
error bars (given in insets only) are calculated as for the MDF 
(see text). 



given by 



Tikpk^ 



+ 0(c^«). (10) 



The simpler scaling for the MDF led us to emphasize this 
quantity in the main body of the text. 



Contrasting the SSF in the DE, GGE and GCE 

In Figs. 8 and 9, we now contrast these results for 
the SSF in the DE with those obtained in the GGE and 
GCE. We plot the results vs. momentum expressed in 
units of kL/kp- The particular form of the SSF in at 
least the diagonal ensemble at small w then suggests that 
in doubling TV and L while halving w, the value of the 



SSF wiU double (see Eqn. [10|). 

For both displayed values of the interaction strength, 
the agreement between the DE and GGE is very good, 
and much better than between the DE and GCE. This 
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FIG. 8: The difference between the SSF as computed in the 
DE and as computed in the GGE and GGE for A^ = L = 28 
(left) and A'^ = L = 56 (right). Top: The gas after release 
from a trap of strength uj — 0.32 (TV — L — 28) and lj = 
0.16 (TV = L = 56) for c = 7200. Bottom: The gas after 
release from the same trap but for c = 10. The GGE estimate 
is seen to be more accurate than the GGE throughout the 
range of values considered. 



is true for the different trap strengths presented in both 
Figs. 8 and 9. We note, however, that the disagreement 
between the DE and GGE is larger for c = 10 than for 
c = 7200 (Fig. 8). 

In terms of a finite size analysis, we note that for the 
data in Fig. 8, the DE-GGE curves maintain, roughly 
speaking, the same shape between the N = L = 28 and 
TV = L = 56. Because of how the SSF is scaling with sys- 
tem size and our choice of units for momenta, this means 
the difference between these two ensembles is decreasing 
as system size grows. However the same cannot be said 
for the DE-GCE curves. For the larger system size, the 
DE-GCE curves are notably more upturned at small mo- 
menta suggesting that in the infinite volume limit, the 
two ensembles will yield different results for the SSF. 

This effect is, however, much less pronounced for the 
data in Fig. 9 where a stronger trap is used. Here the 
DE-GGE curve appears flatter for the larger system size 
data (N=L=56), while, the DE-GCE curve appears much 
the same for the two different system sizes. A more de- 
flnable trend may be elusive here because of the larger 
uncertainties associated with the larger trap values. 



Assessing the convergence of the computation of the 
MDF and SSF in the DE 



MDF 

In computing the MDF in the DE for TV = L = 56 and 
c = 10, we truncated the expression for the ground state 
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And finally for iV = L = 28 and c = 7200, we used 500 
states saturating 0.99997 of the wavefunction norm, with 
f-sum saturations, 0.99998, 0.99998, and 0.99998. 



Calculation of Error Bars 



Diagonal Ensemble 



FIG. 9: The diflterence between the SSF as computed in the 
DE and as computed in the GGE and GCE for TV = L = 28 
(left) and AT = L = 56 (right) for c = 7200, after release 
from traps of strengths iu — 0.512 (A^ — L — 28) and oj = 
0.256 {N = L = 56). 



to 7000 states (for N = L = 14, 100 states were used; 
for iV = L = 28, 500). The sum of the coefficients, |csp, 
over these 7000 states was 0.9884 (for N ^ L ^ U and 
AT = L = 28, respectively 0.9999 and 0.9992). The MDF 
itself can be checked using the sum rule 7; X]fe "fe — T- 
The sum rule saturation for N — L — 56 was 0.9690, for 
AT = L = 14, 0.9995, and for TV = i = 28, 0.9982. 

For TV = L = 56 and c = 7200, we kept 2000 states, 
leading to the truncated sum J^s I'^sP equalling 0.9998. 
The integrated MDF sum rule was saturated to 0.9900. 
For N — L — 14, 100 states were used, yielding satura- 
tions of 0.9999 and 0.9992 respectively. For A^ == L == 28, 
500 states were used, giving 0.9999 and 0.9988. 



SSF 



In computing the SSF in the DE for c = 10 and 
N = L = 56, we again truncated the expression for 
the ground state to the same 7000 states. The f-sum 
rule J^^ ^ujS{k,u!) = j^k^ then provides a convenient 
check on the end result at any fixed momentum. This 
was saturated to 0.9884,0.9766,0.9451 for k = 2tt /L,kF 
and 2kF by using ABACUS for excited states; the lack 
of saturation around 2fci? is visible in the top of Fig. [T] 

For the case of c = 7200 we kept 1000 states in the 
ground state expansion in the basis of Bethe states, yield- 
ing a truncated sum ^^ jcgp equal to 0.9987. This leads 
to a f-sum rule saturation for these three same momenta 
to the same value of 0.9987. 

For the other system sizes, we obtain the following 
saturations. For N — L ^ 14 and c = 10, we used 
100 states saturating 0.9999 of the wavefunction norm, 
and we obtained as the f-sum saturations at the three 
wavevectors, 0.9999, 0.9996, and 0.9991. For iV = L = 14 
and c = 7200, we used 100 states saturating 0.99999 of 
the wavefunction norm, with f-sum saturations, 0.99997, 
0.99997, and 0.99997. For iV = L = 28 and c = 10, 
we used 500 states saturating 0.9992 of the wavefunction 
norm with f-sum saturations, 0.9992, 0.9964, and 0.9915. 



The error bars for the MDF and the SSF in the diago- 
nal ensemble (DE) (as presented in Figs. 1 and 3 of the 
main text, and in Figs. 6-9 of the supplementary ma- 
terial) were obtained according to the following scheme. 
The NRG gives a ground state in the form 



M 



|G^>=E' 



s=l 

where \s) are exact eigenstates of the untrapped gas. This 
sum, formally, should be over all eigenstates in the sys- 
tem (i.e. M = cxd). In running the NRG to determine 
the coefficients c^ , we limit ourselves to a finite number 
of states (approximately M k, 42000 for A^ = L = 56). 
The coefficients arrived at by the NRG always sum to 1, 
i.e. 






= 1. In computing the SSF or MDF in 



the DE, we further truncate this sum in order to make 
the computation of the correlation function numerically 
manageable, for example keeping only the first M' states 
(A/' = 7000 states for A^ = L = 56). The question is 
then what uncertainties these truncations introduce. To 
get a feel for this for the case A^ = L = 56, we compute 
the DE at two different truncation levels, the first using 
all M' — 7000 states. For this truncation the coefficients 
\cs\^ sum to « 0.99 in all cases. But we also look at a sec- 
ond, more severe, truncation where ^ |csp = 0.97. For 
both truncated wavefunctions, we renormalize the wave- 
function so that it has unit norm. We then ascribe the 
uncertainty to our computation of the DE to the differ- 
ence between these two results. 

We readily admit that this is a heuristic but we feel 
that it gives a ball park for the uncertainty and perhaps 
even overestimates it. 

For the cases N = L ^ 14, 28, we feel the convergence 
of the sum ^ jcgp to 1 is so rapid that the results for the 
MDF and the SSF in the DE are completely converged 
and the error due to truncation is negligible (or at least 
smaller than the symbol size used in the plots). 



Generalized Gibbs and Canonical Ensembles 

We first note that in computing the SSF (as well as 
the MDF) in the GGE and GCE, that while the weights, 
g-eo-GGE a^jjfj g-eo-GCE ^ [yi these ensembles are computed 



9 



without fixing the particle number, when we compute the 
trace over states we only perform a partial trace involving 
those states with the same particle number, N. Thus 
the results presented in Figs. 6-9 of the supplementary 
material and 3 of the main text are computed, strictly 
speaking, in fixed-density sub-ensembles. 

To obtain error bars for the GGE and GCE, we com- 
pute the SSF and MDF in two different sub-ensembles 
and take the differences between these computations 
to obtain an (again heuristic) feel for the uncertainty. 
For the N = L = 56 case, the two (sub)-ensembles 
are the same collections of states used for determining 
uncertainties in the DE, one sub-ensemble has a sum 
J2s\'^A'^ ^ 0-99 '^hile the other satisfies J2s\'^s\'^ '^ 0.97. 

For the N ~ L = 28 case, the first sub-ensemble again 
satisfies ^^ jcgp « 0.99 while the second sub-ensemble 
is (differently from the N = L = 56 case), is taken to be 
one half of the states of the first. For this case, keeping 
instead only states that gave a saturation of 0.97 led to 
too few states in the ensemble to compute, even approx- 
imately, the SSF and MDF in the GGE and GCE. 
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